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1.  Introduction 


Since  its  development  almost  20  years  ago,  dissipative  particle  dynamics  (DPD)  ( 1 ,  2)  has 
become  a  robust  tool  for  simulating  soft  condensed  matter  including  polymers,  surfactants,  and 
colloids  (see  e.g.,  Brennan  and  Lisal  [3]  and  references  therein).  Though  often  used  to  simulate 
polymer  systems,  DPD  is  now  being  applied  in  other  areas  such  as  biophysics  ( 4 )  and  shock 
physics  (5).  Moreover,  DPD  models  have  advanced  beyond  the  typical  purely  repulsive  models 
(e.g.,  references  6-9),  including  augmentations  that  capture  noncentral  shear  forces  (7,  10). 

When  applying  the  DPD  method,  numerical  integration  of  the  equations-of-motion  (EOM)  is  a 
key  consideration,  where  the  stochastic  component  requires  special  attention.  Efficient  and 
accurate  integration  schemes  are  required  to  maintain  reasonable  time  step  sizes,  thus  allowing 
for  the  simulation  of  mesoscale  events.  Moreover,  the  advent  of  conservative  forces  extending 
beyond  purely-repulsive  models  that  contain  attractive  character  further  supports  the  need  for 
effective  integration  schemes.  However,  the  integration  is  a  nontrivial  task  due  to  the  pairwise 
coupling  of  particles  through  the  random  and  dissipative  forces  (11).  Furthennore,  self-consistent 
solutions  are  often  necessary  because  the  dissipative  forces  and  the  relative  velocities  of  the 
particles  are  interdependent,  where  the  considerable  computational  cost  associated  with  this  task 
has  motivated  the  development  of  various  integration  schemes.  Recent  independent  studies  by 
Nikunen  et  al.  (11)  and  Chaudhri  and  Lukes  (12)  assessed  the  quality  and  performance  of  several 
applicable  integration  schemes  for  constant-temperature  DPD,  specifically,  the  velocity- Verlet- 
based  integrator  (13),  self-consistent  Pagonabarraga-Hagen-Frenkel  integrator  (14),  self- 
consistent  velocity- Verlet  integrator  (15,  16),  den  Otter-Clarke  integrator  (17),  first-  and  second- 
order  Shardlow-splitting  algorithms  (SSAs)  (18),  and  Lowe-Andersen  integrator  (19).  Still  other 
schemes  are  possible  such  as  the  Peters  integrator  (20)  and  a  scheme  similar  to  the  SSA 
developed  by  Litvinov  et  al.  (21)  for  highly  dissipative  smoothed-particle  dynamics.  Further, 

De  Fabritiis  et  al.  (22)  developed  numerical  integrators  for  stochastic  models  using  the  Trotter 
expansion,  and  later  applied  the  stochastic  Trotter  expansion  to  the  constant-temperature  DPD 
method  (23).  They  tested  the  DPD-Trotter  integration  scheme  on  a  DPD  fluid  and  found 
accuracy  comparable  to  the  SSA  but  with  a  10%  speed-up.  The  consensus  from  all  of  these 
studies  identified  the  SSA  (18)  as  the  best-performing  approach. 

The  SSA  decomposes  the  EOM  into  differential  equations  that  correspond  to  the  deterministic 
dynamics  and  the  elementary  stochastic  differential  equations  that  correspond  to  the  stochastic 
dynamics.  In  the  original  SSA  formulation,  both  types  of  differential  equations  are  integrated  via 
the  velocity- Verlet  algorithm  (18),  where  the  stochastic  dynamics  are  additionally  solved  in  an 
implicit  manner  that  conserves  linear  momentum.  Previously,  the  SSA  had  been  applied  only  to 
the  constant-temperature  DPD  method,  while  an  SSA-inspired  approach  was  implemented  for  a 
method  in  the  spirit  of  the  energy-conserving  DPD  method  (5). 
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In  this  work,  we  modify  the  original  SSA  formulation  for  constant-temperature  DPD  systems  of 
unequal-mass  particles,  followed  by  the  development  of  an  alternative  algorithm  to  implicitly 
integrate  the  stochastic  dynamics.  The  SSA  is  formulated  for  both  a  velocity- Verlet  scheme  and 
an  implicit  scheme,  where  the  velocity-Verlet  scheme  consistently  performed  better.  For 
completeness,  the  derivations  of  the  Fokker-Planck  equation  (FPE)  and  the  fluctuation- 
dissipation  theorem  (FDT)  are  included. 


2.  Formulation  of  DPD  at  Fixed  Temperature  Using  Shardlow-Like  Splitting 
Numerical  Discretization 


2.1  General  Formulation  of  DPD 

DPD  particles  are  defined  by  a  mass  mi ,  position  r; ,  and  momentum  p, .  The  particles  interact 
with  each  other  via  a  pairwise  force  F/y  that  is  written  as  the  sum  of  a  conservative  force  F,'  , 
dissipative  force  F" ,  and  random  force  F*  : 

i  If  '  'V  0) 

F(/r  is  given  as  the  negative  derivative  of  a  coarse-grain  potential  u]-  ' ,  expressed  as 


r  d  u  a 
pt  _  _ lJ 
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CG 


r. 
v  y 


(2) 


where  ry  =  r;.  -  r;  is  the  separation  vector  between  particle  i  and  particle  j ,  and  /y  =  |r;;  | .  The 

remaining  two  forces,  Ff  and  Ff  ,  can  be  interpreted  as  a  means  to  compensate  for  the  degrees 

of  freedom  neglected  by  coarse  graining.  The  conservative  force  is  not  specified  by  the  DPD 
formulation  and  can  be  chosen  to  include  any  forces  that  are  appropriate  for  a  given  application, 
including  multibody  interactions  (e.g.,  references  6,  8,  and  24]).  and  F*  are  defined  as 


and 


F/f  =  aijC0 


!hk 


rjL 

ru 


(3) 

(4) 


where  ytj  and  afj  are  the  friction  coefficient  and  noise  amplitude  between  particle  i  and 

p  P 

particle  j  ,  respectively,  and  v;  =  — - —  and  W..  are  independent  Wiener  processes  such  that 

m,  in 

Wj  =  Wjf .  The  weight  functions  coD(r)  and  o/:(r)  vanish  for  r>rc,  where  rc  is  the  cutoff  radius. 
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Note  that  F5  is  completely  independent  of  F(J  and  Ff  ,  while  F/J  and  F*  are  not  independent 

but  rather  are  coupled  through  a  fluctuation-dissipation  relation.  This  coupling  arises  from  the 
requirement  that  in  the  thermodynamic  limit,  the  system  samples  the  corresponding  probability 
distribution.  The  necessary  conditions  can  be  derived  using  an  FPE  (25,  26),  which  is  presented 
in  this  work. 


2.2  Constant-Temperature  DPD 

For  constant-temperature  and  constant-volume  conditions,  the  evolution  of  DPD  particles  in  time 
t  is  governed  by  the  following  EOM: 


dr(.  =  —  dt 
m . 


(5a) 


dp-=Z(F^+F^+F"^ 


(5b) 


where  N  is  the  total  number  of  particles.  This  EOM  set  conserves  total  momentum  P  -  ^  p  as  a 

i 

consequence  of  pairwise  additivity  of  the  conservative,  dissipative,  and  random  forces. 

Espanol  and  Warren  (27)  showed  that  this  system  samples  the  canonical  probability  distribution 
and  obeys  the  FDT  if  the  following  relations  hold: 


c o 


=  2r,jkJ 


(6) 


where  kB  is  the  Boltzmann  constant  and  T  is  the  temperature.  For  completeness,  the  FPE  and  an 
outline  of  the  derivation  of  the  FDT  is  provided  in  appendix  A.  coD(r)  and  coR(r)  are  typically 
chosen  (13)  as 


con(r)=  coR(r)~ 


( 

1 

v 


=  0 


(r  <  r  ) 
[r>rc) 


(7) 


Note  that  the  dependence  of  the  random-force  noise  amplitude  cr.  on  the  temperature  T  can  be 
physically  interpreted  as  the  system  being  in  thermal  contact  with  a  heat  reservoir. 


2.2.1  Numerical  Discretization 

Next  we  review  the  first-order  SSA  (18).  The  original  SSA  was  formulated  for  systems  of 
particles  with  equal  masses  (18),  while  here  we  provide  a  fonnulation  for  systems  with  unequal 
masses.  The  basic  idea  of  the  SSA  is  the  decomposition  of  the  EOM  (equations  5a  and  5b)  into 
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differential  equations  corresponding  to  the  detenninistic  dynamics  and  elementary  stochastic 
differential  equations  (SDEs)  corresponding  to  the  fluctuation-dissipation  contributions.  In  the 
original  SSA  formulation,  integration  of  both  types  of  differential  equations  uses  the  velocity- 
Verlet  algorithm  (28). 

For  the  constant-temperature  DPD  SSA  fonnulation,  the  conservative  tenns  consist  of  the 
following  differential  equations: 

dr(  =  —  dt  (8a) 

mi 

(i  =  l . N)  , 

dp,=2Xd(  (8b) 

j*i 


while  the  fluctuation-dissipation  tenns  consist  of  the  elementary  SDEs: 

f  V. 


dp':'  =  -rM 


dp  =  dp 
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r.. 
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for  each  i  <  j 


(9) 


where  the  superscript  i  -  j  indicates  that  the  variation  of  momenta  is  considered  for  a  pair  of 
interacting  particles  i  and  j  only,  and  dWtj  =  dWjt  are  the  increments  of  the  Wiener  processes. 
Note  that  equation  5b  is  recovered  by  adding  ^dpp7  to  equation  8b. 

j*i 


Both  the  conservative  and  the  fluctuation-dissipation  terms  can  be  solved  using  velocity- Verlet 
algorithms,  where  the  SSA  occurs  in  the  following  manner.  First,  denote  the  overall  solution 
operator  <j)M  as  the  stochastic  flow  map  for  a  time  step  At  corresponding  to  the  dynamics 

defined  in  equations  5a  and  b,  and  further  denote  (/)(M  and  .  (l  <  /  <  j  <  N)  as  the  solution 

operators  or  flow  maps  associated  with  equations  8a  and  b,  and  equation  9,  respectively. 
(Stochastic  flow  maps  are  mathematical  constructs  that  define  solutions  to  SDEs,  where  the 
solution  operators  are  performed  in  a  successive  manner.  See  Shardlow  (18)  and  references 
therein  for  further  discussion  of  stochastic  flow  maps  and  splitting  methods  for  SDEs.)  An 
approximation  of  (f)y  for  first-order  splitting  (18)  can  then  be  given  as 
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where  °  denotes  a  composition  of  operators. 

For  each  fluctuation-dissipation  term  ,  momenta  can  be  updated  using  the  velocity- Verlet 
scheme  (18): 
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where  the  superscript  i  -  j  has  been  omitted  for  notational  simplicity,  and  gf.  =  gjt  is  a  Gaussian 
random  number  with  zero  mean  and  unit  variance  that  is  chosen  independently  for  each  pair  of 
interacting  particles.  (The  use  of  Va t  rather  than  At ,  in  equations  1  la-d,  is  the  consequence  of 
a  discrete  Wiener  process.)  Equations  lie  and  d  are  implicit  due  to  the  presence  of  the 
v(/(t  +  At)  tenns  on  the  right  side  of  each  equation.  However,  they  can  be  rewritten  in  an  explicit 

form  (18)  that  eliminates  these  tenns,  which  are  then  given  as: 
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As  an  alternative  to  the  SSA  fonnulation  that  detennines  the  fluctuation-dissipation  tenns  using 
a  velocity- Verlet  algorithm,  here  we  also  provide  a  modified  SSA,  where  the  velocity- Verlet 
algorithm  is  replaced  by  an  implicit  algorithm  based  upon  the  following  arguments.  For 
conservative  systems,  the  Liouville  equation  indicates  that  the  flow  of  system  points  in  phase- 
space  is  incompressible.  The  velocity- Verlet  algorithm  is  reversible  and  area-preserving,  and 
thus  possesses  this  desirable  property  for  molecular  dynamics  simulations  (29-32).  However, 
DPD  incorporates  dissipative  and  random  forces  such  that  the  volume  element  in  phase-space  is 
not  preserved  and  the  trajectories  are  not  time-reversible.  Moreover,  the  motion  of  the  density  of 
system  points  in  phase-space  is  a  Fokker-Planck  (diffusion)  equation,  leading  to  a  physical 
increase  of  the  phase-space  volume  as  the  diffusive  dynamics  irreversibly  evolves.  Therefore, 
ensuring  time  reversibility  of  the  integration  of  is  not  necessary.  Consequently,  we  propose 

the  following  implicit  algorithm  as  an  alternative  to  the  velocity- Verlet  algorithm  to  integrate 
this  contribution: 


Analogous  to  the  derivation  of  equations  lie  and  Ilf,  equations  1  lg  and  1  lh  can  be  rewritten  in 
an  explicit  fonn  as: 
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Finally,  the  conservative  term  <f>^  is  approximated  using  the  velocity-Verlet  algorithm,  where  the 
momenta  and  positions  are  updated  as 
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p,.(t  +  A?)  =  p, 


r  At  At 
t  H - 

V  2) 


+  — F,c(f  +  Af)  (i  =  l,...,N) 


(12) 


We  denote  the  original  SSA  given  by  equation  10,  equations  1  la,  1  lb,  1  le,  1  If,  and  12  as 
SSA-VV,  while  the  modified  SSA  given  by  equation  10,  equations  1  li  and  1  lj,  and  equation  12 
as  SSA-I.  The  following  outlines  summarize  the  SSA-VV  and  SSA-I  in  a  practical  form.  By 
comparing  these  forms,  it  is  evident  that  the  SSA-I  requires  fewer  mathematical  operations  than 
the  SSA-VV  and  thus  is  less  computationally  intensive.  However,  as  demonstrated  in  section  IV, 
the  SSA-VV  is  a  better  performing  integration  scheme.  Finally,  note  that  for  the  proposed 
splitting  schemes  in  the  limit  of  vanishing  dissipative  and  random  forces,  the  reversibility  of  the 
trajectories  and  the  area-preserving  property  is  recovered,  as  expected. 


(I)  Outline  of  the  SSA-VV  for  constant-temperature  DPD: 

1 .  Stochastic  Integration  for  all  i  -  j  pairs  of  particles: 
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2.  Deterministic  Integration  no.  1  for  i  =  1  : 
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3.  Conservative  Force  Calculation,'.  |fc  };3, 

4.  Deterministic  Integration  no.  2  for  i  =  1,..., N : 

A  t  ~ 


(II)  Outline  of  the  SSA-I  for  constant-temperature  DPD: 

1 .  Stochastic  Integration  for  all  i  -  j  pairs  of  particles: 
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2.  Deterministic  Integration  no.  1  for  i  =  \,...,N : 
_  .  At  c 


(a)  P, 

(b)  r, 


■  r .  +  At  — 


m 


{c 

F(.  ).=1 

4.  Deterministic  Integration  no.  2  for  i  =  1,..., N : 
_  ,  _  .  At  c 


3.  Computational  Details 


The  SSA-VV  and  SSA-I  for  the  constant-temperature  DPD  method  were  tested  using  the 
standard  DPD  fluid  (13),  where  details  of  the  conservative  forces  for  this  model  are  given  in 
appendix  B.  Both  a  pure  component  case  and  an  equimolar  binary  mixture  were  tested  for  the 
DPD  fluid  model,  with  system  sizes  of  Af  =  10125.For  these  simulations,  the  following  reduced 
units  were  used:  rc  is  the  unit  of  length;  the  mass  of  a  DPD  particle  is  the  unit  of  mass;  and  the 

unit  of  energy  is  kBT  .  Using  these  reduced  units,  we  set  the  maximum  repulsion  between 
particles  i  and  j  as  at]  =  25  for  the  pure  DPD  fluid,  and  as  <x  =  25  and  28  for  the  like  and 
unlike  i  -  j  interactions,  respectively,  for  the  binary  DPD  fluid.  Further,  for  all  cases,  we  set  the 
noise  amplitude  cr.  =  3 . 
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4.  Results  Comparison  of  SSA-VV  and  SSA-I  Performance 


For  both  the  pure  and  equimolar  binary  DPD  fluids,  we  performed  constant-temperature  DPD 

N 

simulations  at  particle  density  p  =  —  =  3  and  T  =  1  using  the  SSA-VV  and  SSA-I,  where  At 

was  varied  from  0.01  to  0.15.  For  the  various  time  steps  and  runs  of  length  trun  =  5000,  figure  1 
(pure  fluid)  and  figure  2  (equimolar  binary  fluid)  compare  values  of  the  simulated  kinetic 
/  1  "  n..n.\ 

temperature  (Tkin)  =  ( - V  — '■ — )  to  the  prescribed  T ,  along  with  the  configurational  energy 

\3 kBN  <=1  mt  / 

per  particle  (u) ,  and  virial  pressure  (Pvir)  ((.)  denotes  an  ensemble  average).  Values  of 
{{Tkin)-  t)/T  ,  (u)  and  (Pvir)  for  At  >0.07  begin  to  significantly  diverge  and  therefore  are  not 

shown  in  figures  1  and  2.  For  both  fluids,  the  SSA-I  systematically  underestimates  T  by 
approximately  l%-2%  for  At  below  0.07,  but  more  importantly  its  performance  shows  non¬ 
monotonic  dependence  on  At ,  which  is  undesirable  behavior  for  a  numerical  integrator.  From  a 
practical  viewpoint,  figures  1  and  2,  and  the  values  of  other  properties  such  as  the  radial 
distribution  function  and  self-diffusion  coefficient  (not  shown  here),  suggest  that  for  all  values  of 
At  <  0.05  ,  the  results  of  the  constant-temperature  DPD  simulations  are  within  good  agreement 
for  both  the  SSA-W  and  SSA-I. 

For  completeness,  also  shown  in  figures  1  and  2  is  the  standard  velocity- Verlet  (SVV)  algorithm 

(13)  often  used  for  DPD  simulations.  Comparing  the  SSA-VV  and  SVV,  we  see  that 

(( Tki„ )~  t)/T  is  positive  and  increases  with  increasing  At  for  both  algorithms.  However,  the 

SSA-VV  is  superior  in  comparison  with  the  SVV  for  both  the  pure  (in  agreement  with  previous 
work  [  /  A] )  and  binary  DPD  fluids. 
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Figure  1.  (a)  The  simulated  kinetic  temperature  y  kin  )  and  (b)  the  configurational  energy 

per  particle  (w)  and  the  virial  pressure  ( Pvir }  of  the  pure  DPD  fluid  as  a  funct;on 
of  the  integration  time  step  At  for  constant-temperature  DPD  simulations  with  the 
SSA-VV,  SSA-1,  and  SVV  integrators  at  p  —  3  and  T  =  1 
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(a) 


(b) 
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Figure  2.  (a)  The  simulated  kinetic  temperature  (Tkin^  and  (b)  the  configurational  energy  per 

particle  (ll^  and  the  virial  pressure  of  the  equimolar  binary  DPD  fluid  as  a 

function  of  the  integration  time  step  fa  for  constant-temperature  DPD  simulations 
with  the  SSA-VV,  SSA-I,  and  SVV  integrators  at  p  =  3  and  T  =1 . 
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5.  Conclusion 


A  comprehensive  description  of  numerical  integration  schemes  based  upon  the  Shardlow- 
splitting  algorithm  was  presented  for  the  constant-temperature  DPD  method.  The  original  SSA 
formulated  for  systems  of  equal-mass  particles  has  been  extended  to  systems  of  unequal-mass 
particles.  Both  a  velocity- Verlet  scheme  and  an  implicit  scheme  were  formulated  for  the 
integration  of  the  fluctuation-dissipation  contribution,  where  the  velocity- Verlet  scheme 
consistently  performed  better. 
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The  Fokker-Planck  equation  (FPE)  corresponding  to  the  equations  of  motion  given  by  equations 
5a  and  5b  of  this  report  is 

~T~  =  Lcp  +  LDp,  (A-l) 

dt 


where  p  =  /?(r,p;t)  is  the  time-dependent  probability  density,  r  =  {r.}^  and  p  =  {pi  are 
particle  positions  and  momenta,  respectively,  and  N  is  the  total  number  of  particles.  In  equation 
A- 1 ,  we  have  defined  the  operators 


Lc  =~l 


Z-—+ZZF, 

V  /  rnt  or,  ,,,  ■ 


5p, 
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**-ZZ 
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(A-3) 


where  mi  is  the  mass  of  particle  i ,  r,y  =  r;.  -  r  is  the  separation  vector  between  particle  i  and 
particle  j  ,  rtj  =  r(/ 1 ,  F,y  is  the  conservative  force  acting  between  particle  i  and  particle  j  ,  ytj 
and  cr.  are  the  friction  coefficient  and  noise  amplitude  between  particle  i  and  particle  j , 


respectively,  v,  =  — -  —  ,  and  coD  and  coR  are  weight  functions  of  the  dissipative  and  random 
mi  in 

forces,  respectively.  Lc  is  the  Liouville  operator  corresponding  to  a  Hamiltonian  system 
interacting  with  conservative  forces  F  .  The  operator  LD  takes  into  account  the  effects  of  the 
dissipative  and  random  forces. 


The  equilibrium  probability  density  peq  =  pe?(r,p)  should  be  a  steady-state  solution  of  the  FPE, 
dp 

—  =  0  .  Thus,  the  FDT  can  be  derived  by  imposing  that  peq  corresponds  to  the  canonical 

dt 

(Gibbs-Boltzmann)  probability  density: 


Peq{*,  p)  =  -^-exp 


KT 


=  zexp 
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kj 


(A-4) 
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where  //(r,p  )  is  the  Hamiltonian  of  a  dissipative  particle  dynamics  (DPD)  system,  Z  is  the 
nonnalizing  partition  function,  kB  is  the  Boltzmann  constant,  T  is  the  system  temperature,  and 
UyG is  a  coarse-grain  potential.  Lcpeq  =  0  since  the  canonical  probability  density  is  the 

equilibrium  solution  for  the  conservative  system.  The  FDT  then  follows  from  the  requirement 
LDpeq  =  0  ,  which  is  satisfied  for 


i.e.,  for 


=  iKjhT 
(a>R  y = cod 


(A-5) 


The  physical  interpretation  of  the  FDT,  equation  A-5,  is  the  following.  In  a  constant-temperature 
DPD  simulation,  the  system  is  implicitly  connected  to  a  heat  reservoir  whose  temperature 
corresponds  to  the  imposed  temperature.  During  the  simulation,  heat  is  exchanged  between  the 
heat  reservoir  and  the  system  such  that  the  average  temperature  of  the  system  equals  the 
reservoir  temperature.  Thus,  the  FDT  indicates  that  the  amplitude  of  the  random-force  friction 
coefficient  <x.  depends  upon  the  system  temperature. 
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For  the  models  considered  in  this  work,  the  details  of  the  conservative  forces  expressed  in 
equation  2  of  the  main  text  are  the  following:  uGG  for  the  pure  and  binary  dissipative  particle 

dynamics  fluids 1  is  given  by 

«r=vy>M  (B-i) 

where  is  the  maximum  repulsion  between  particle  i  and  particle  j  . 


'Groot,  R.  D.;  Warren,  P.  B.  Dissipative  Particle  Dynamics:  Bridging  the  Gap  Between  Atomistic  and  Mesoscopic 
Simulation.  J.  Chem.  Phys.  1997,  107,  4423. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


DPD 

dissipative  particle  dynamics 

EOM 

equations-of-motion 

FPE 

Fokker-Planck  equation 

FDT 

fluctuation-dissipation  theorem 

SDEs 

stochastic  differential  equations 

SSA 

Shardlow-splitting  algorithm 

SSA-VV 

Shardlow-splitting  algorithm-velocity  Verlet 

SSA-I 

Shardlow-splitting  algorithm-implicit 

SVY 

standard  velocity- Verlet 
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